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Abstract. We have performed a statistical analysis using broad 
band IRAS data on about 500 planetary nebulae with the aim 
of characterizing their dust content. Our approach is different 
from previous studies in that it uses an extensive grid of pho- 
toionization models to test the methods for deriving the dust 
temperature, the dust-to-gas mass ratio and the average grain 
size. In addition, we use only distance independent diagrams. 

With our models, we show the effect of contamination by 
atomic lines in the broad band IRAS fluxes during planetary 
nebula evolution. We find that planetary nebulae with very dif- 
ferent dust-to-gas mass ratios exist, so that the dust content is a 
primordial parameter for the interpretation of far infrared data 
of planetary nebulae. In contrast with previous studies, we find 
no evidence for a decrease in the dust-to-gas mass ratio as the 
planetary nebulae evolve. We also show t hat th e decrease in 
grain size a dvocated by Natta & Panagia ( 1981 ) and Lenzuni 
et al. (1989) is an artefact of their method of analysis. Our re- 
sults suggest that the timescale for destruction of dust grains in 
planetary nebulae is larger than their lifetime. [] 

Key words: planetary nebulae: general - dust - Stars: AGB 
and post- AGB - circumstellar matter 



1. Introduction 

Far infrared measurements of planetary nebulae allow one to 
investigate their dust content. From broad band photometry, 
one can derive the characteristic temperature of the emitting 
dust and the total amount of energy that it radiates. 

For example, Natta & Pan agia ( 198 1 ), using the aircraft ob- 
servations by Moseley ( 1980 ), analyzed a sample of 10 plane- 
tary nebulae, and found evidence that the dust-to-gas mass ra- 
tio, md/wg, and the grain size were decreasing systematically 
with nebular radius, indicating a time variation of t he dust prop- 
erties in planetary nebulae. Then, Pottasch et al. ( 1984a| ) used 
IRAS broad band photometry to analyze a sample of 46 plane- 
tary nebulae. They also found rrid/mg to decrease with radius. 

Send ojfprint requests to: Grazyna Stasinska (grazyna@obspm.fr) 
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being as large as a few 10 ^ for the smallest nebulae and two 
order of magnitudes smaller for the largest ones. Later, Lenzuni 



et al. (1989) repeated the analysis of Natta & Panagia (1981) 



on a much larger sample: they analyzed 233 planetary nebu- 
lae which had been observed by IRAS in at least the 25 jim. 
and 60 /im bands and appeared in the Daub ( 1982 ) catalog of 
distances . Esse ntially, they confirmed the finding of Natta & 
Panagia ( |l98l| ). 

The problem of the infrared luminosity of planetary neb- 
ulae has been addressed by several authors: Pottasch et 
al. (|l984a| ), Iyengar (|1986| ), Zijlstra et al. (|1989| ), Ratag et 
al. ( |199C| ), Amnuel ( 11994 . Pottasch et al. ( |l984ij ) found that 
very high infrared excesses (IRE, i.e. the infrared luminosity in 
terms of what can be accounted for by absorption of the Lya ra- 
diation produced in the nebula) were usually found in the small, 
high dust temperature nebulae. They concluded that, in these 
cases, dust is probably mostly heated by the radiation of the 
central star on the longwave side of Lya, since the nebulae are 



young. Ratag et al. (1990) also considered the possibility that 



high IRE could be due to high dust contents or to heating by 



the interstellar radiation field. Amnuel (1994) favored the latter 
explanation. 

All the studies above based their conclusions either on a 
simple, empirical interpretation of the data, or, in the best cases, 
on semi-analytical models of planetary nebulae that involve 
drastic simplifications. Also, excepting the studies specifically 
devoted to planetary nebulae in the Galactic bulge (Ratag et 
al. 1990), they rely heavily on the adopted distances to the neb- 
ulae. Therefore, they must be considered with caution. A safer 
approach, in our opinion, for studies of statistical nature, is to 
interpret the observed data with the help of grids of numeri- 
cal models that are able to compute the observables and use 
distance independent diagrams. Such an approach for the inter- 
pretation of the infrared emission of planetary nebulae has been 



attempted by Volk (1992i and Vil'koviskii & Efimov (1992) 



Unfortunately, both these studies were rather exploratory, and 
did not address explicitly the questions mentioned above. 

Fifteen years after the publication of the results of the IRAS 
mission, there is still work to be done on their interpretation 
as regards the statistical dust properties of planetary nebulae. 
It is the purpose of the present paper to contribute to this, by 
analyzing the broad band IRAS fluxes of planetary nebulae in 
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the light of series of dusty photoionization models for evolv- 
ing planetary nebulae using distance independent diagrams. In 
Section 2, we briefly present the observational material. In Sec- 
tion 3, we describe the photoionization models. In Section 4, 
we analyze the contamination of IRAS broad band fluxes by 
line emission. In Section 5, we discuss the dust temperature. 
Section 6 is devoted to the dust-to-gas mass ratio, Section 7 to 
the infrared excess and Section 8 to the grain size. A general 
discussion of the main findings is given in Section 9. 



2. The observational data base 

We have considered the planetary nebulae from the Strasbourg- 
ES O Ca talogue of Galactic Planetary Nebulae (Acker et 
al 



restrictive conditions on the IRAS fluxes and more recent com- 



pilation of H/3 fluxes and angular diameters than Daub (1982) 



1992| ) that appear in the IRAS Point Source Catalog ( |1988D 
and have flux quality Q = 3 in at least the 25 and 60 /im bands. 
According to the IRAS Point Source Catalog, the error in the 
fluxes with (3 = 3 is about 10%. There are 548 such objects. 
The inclusion of planetary nebulae with Q = 2 adds only 86 
objects with lower quality data, and does not change the con- 
clusions of this paper 

Since individually determined distances in planetary neb- 
ulae are scarce and are sometimes contradictory for the same 
object when using different methods (e.g. Sabbadin 1986), we 
consider only distance-independent diagrams. We use the red- 
dening corrected H/3 fluxes (or, equivalently, the radio fluxes at 
6 cm converted into H/3) and angular diameters from the same 
sources as Tylenda & Stasinska (1994 ). We discard data that are 
only upper limits. We end up with 477 nebulae for which we 
can compute the average surface brightness, ^(H/?), and which 
have the required IRAS data (i.e. 25 and 60 ^m fluxes with 
Q = 3). Table 1^ lists, in columns 3 and 4, the adopted values 
of the reddening corrected H/3 fluxes and angular radii for all 
the nebulae having 25 and 60 fj,m fluxes with Q = 3 (columns 1 
and 2 give the planetary nebula PN G identifications and usual 
names, respectively). As is known, and seen from our models 
below, S(li(3) is a rather robust estimator of the planetary neb- 
ula age. It is safer than the estimated planetary nebula radius, 
which was used in former studies on the dust-to-gas mass ratio 
in planetary nebulae, among other reasons, because it is dis- 
tance independent. 

Finally, to discuss the dust-to-gas mass ratio, we need an es- 
timate of the gas density. We take the electron densities derived 
from the [S II] ratios as compiled by Gorny (private communi- 
cation), again discarding upper or lower limits. These densities 
are listed in Column 5 of Table 1 . The main references are the 



compilation of Stanghellini & Kaler (1989), and the observa- 
tions reported by Acke r et al . ( 1989 ), Acker et al. (1991) and 
Kingsburgh & Barlow (1994). Thus, we end up with a sample 
of 230 objects for which we can discuss the time evolution of 
md/ mg. Th is sample is similar in size to the one of Lenzuni 



et al. (1989), but differs slightly from it because of our more 



Table 1 is accessible in electronic form at the CDS via anony- 
mous ftp to cdsarc.u-st rasbg.fr (130.79.128.5) or via h ttp://cdsweb.u- 
strasbg.fr/Abstract.html 



3. The photoionization models 

The modelling follows the same policy as in Stasinska et 
al. (1998). We construct sequences of photoionization mod- 



els that represent the evolution of a nebula around a central 
star that has left the AGB and is progressing towards the white 
dwarf stage. The nebula is modelled as a uniformly expand- 
ing spherical shell of given total mass (with a relative thickness 
(^^out — Rin) /Rout = 0.3, where i?in and i?out are the inner and 
outer radius of the nebula, respectively) in which the density is 
supposed uniform. In such a simple representation, the param- 
eters defining a planetary nebula are the central star mass, Af*, 
the nebular mass, M^eb^ its expansion velocity, Vcxp, and its 
chemical composition (also assumed uniform). True nebulae 
are, of course, more complicated than such a simple descrip- 
tion. More elaborate dynamical models have been computed, 
even includi ng de parture from spherical symmetry (e.g. Frank 
& Mellema 1994 ), which aim at a better representation of the 
real evolution of planetary nebulae (but those models do not 
include dust and only a few sequences are computed). Our ap- 
proach has the virtue of being sufficiently simple to allow one 
to run series of models covering the whole parameter space 
suggested by the observations and to make easier the interpre- 
tation of observational data, in the framework of general ideas 
on the evolution of planetary nebulae. 

The models are computed with the photoionization code 
PHOT O, usin g the atomic data described in Stasinska & Lei- 
therer (1996). The code assumes spherical geometry, with a 
central ionizing source. The diffuse radiation is treated assum- 
ing that all the photons are emitted outwards in a solid angle of 
2tt, and the transfer of the resonant photons of hydrogen and 
helium is computed with the same outward only approxima- 
tion, but multiplying the photo-absorption cross-section by an 
appropriate factor to account for the increased path lenght due 
to scattering (Adams 1975). In contrast with the model plane- 
tary nebulae presented in Stasinska et al. (1998), which were 
aimed at discussing the emission from the ionized region only, 
the present models must take into account the possible pres- 
ence of a neutral zone surrounding the ionized region, where 
dust emission also occurs. This is done by extending the com- 
putation of the radiation transfer outside the ionized zone, until 
the total nebular mass assigned to a given model is reached. The 
gas temperature is not computed in the neutral zone, but simply 
set to lOOK, since it has no incidence on the results presented 
here. 



Dust is introduced as described in Harrington et al. (1988) 



In short, it is assumed to be composed of spherical grains with 
the classical size distribution (oc a~^^^, where a is the grain 



size) of Mathis et al. ( 1977 ). At each nebular radius, the equa- 



tion of radiative equilibrium is solved, for 18 grain sizes log- 
arithmically spaced between 0.005 and 0.25 fj,m. Dust in the 
models is assumed to interact only with the photons, not with 
the gas particles. It is heated by the attenuated stellar radiation 
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field and the diffuse radiation field produced in the nebula. At 
each radius, once the grain temperature is known for each grain 
size, the infrared emissivity is computed at 256 wavelenghts. 
The total luminosity at each wavelength is calculated by inte- 
gration over the nebular volume. The luminosities in the IRAS 
bands are obtained by convolution with the response function 
of the IRAS filters. As a check, the same procedure has been 
applied to the photoionization code written by Szczerba (see 



e.g. Ggsicki et al. 1996). Both codes give consistent results. 



and are ab le to r eproduce the models of IC418 and NGC7662 
by Hoare ( |l990| ). 

The central stars are assumed to evolve according to the 
H-burning post-AGB evolutionary tracks as interpolated by 



Gorny et al. ( 1 997 ) from the tracks of Blocker ( 1 995 ) and refer 



ences therein. The stars are assumed to radiate as blackbodies. 

Results are presented with only one nebular abun- 
dance set (He/H = 0.084; C/H = 8.210-5; N/H = 2.210-^ 
O/H = 4.2 10-4; Ne/H = 8.5 IQ-^; Mg/H = 2. 1 IQ-^; 
Si/H = 2.2 10-5; S/H= 1.0 10-5;C1/H= 1.6 10^^; 
Ar/H=3.010-^;Fe/H = 2.410-S), corresponding to half 
solar metallicity and depleted metals. In the present study, 
the only important effect of the gas chemical composition is 
to affect the energy emitted in the lines included in the IRAS 
bands (see below). 

IRAS Low Resolution Spectroscopy of planetary nebu- 
lae (see e.g. Pottasch 1987) has shown that planetary nebu- 
lae contain several kinds of dust grains. Generally, it is be- 
lieved that, in carbon rich planetary nebulae, dust is com- 
posed of carbon-based grains, e.g. graphite, anthracite, amor- 
phous carbon and/or policyclic aromatic hydrocarbons, while 
in oxygen-rich planetary nebulae it is composed of oxygen- 
based grains, e.g. different forms of silicates. However, ob- 
servations with the Infrared Space Observatory showed that 
some nebulae around Wolf-Rayet central stars have both fea- 
tures (carbon- and oxygen-based) in their spectra (Waters et 
al. 



1998) 



The optical properties of the various materials are differ- 
ent and this affects the infrared emission from the dust. The 
number of planetary nebulae in which the carbon abundance 



has been estimated is only about 50 (Rola & Stasinska 1994) 



and even in these objects, it is often very uncertain. Therefore, 
in our statistical approach, we must take into account the fact 
that the grain composition varies from object to object in an 
unknown manner We have therefore computed models with 



different grain compositions: graphite (Drai ne & Laor 1993), 



amorphous carbon (Rouleau & M artin |1991[ ) and circumstellar 
silicates (David & Pegourie 1995 ). 

In view of the large parameter space to consider, we have 
chosen a reference model with A/* = 0.60 M©, Mnch = 0.2 M0, 
t^exp = 20kms^^ and dust grains composed of graphite, with 
md/?Tig = 7.5 10-4 (which corresponds to a dust-to-hydrogen 
mass ratio TOd/"^H = 10^"^). In the figures presented below, 
we discuss differences with respect to this reference model by 
changing only one defining parameter. In the panels with la- 
bel a, takes the values 0.58, 0.60, 0.62 and 0.64 Mg; in 
panels b. A/nob = 0.1, 0.2 and 0.4 M©; in panels c, fcxp= 10, 



20 and 40kms-i; in panels d, TOd/mH = 10-4, 10-^ 10-^ 
and 10^4; finally, in panels e we explore various dust compo- 
sitions: graphite, amorphous carbon and circumstellar silicates. 
It is then relatively easy, from these figures, to imagine the be- 
havior of a model of any arbitrary combination of defining pa- 
rameters. 

All the sequences of models have been computed with a 
time step of 500 yr, until an age of lOOOOyr was reached, and 
each symbol in the figures below represents one model. There- 
fore zones where points accumulate in a theoretical diagram 
correspond to zones where observational points should accu- 
mulate too (if observational selection does not discriminate 
against this particular zone). It would have been more satis- 
factory to build the models until a given surface brightness was 
reached. This, however, would result in very long sequences 

for some cases, since it is easily shown that the time to reach a 

2/5 _ 1 

given surface brightness goes like M^^^v^^^ when the nebulae 
have reached the optically thin regime. 

4. The contamination of IRAS broad band fluxes by line 
emission in the model nebulae 



As discussed in Pottasch et al. ( 19844 b) and Lenzuni et 
al. (1989), the broad band fluxes measured by IRAS are not 



only due to dust emission. There is a contribution from atomic 
line emission, as seen directly in the low resolution IRAS spec- 
tra of planetary nebulae (see e.g. Pottasch et al. 1984b). Be- 



cause this effe ct was e specially important in the 12/im band, 
Pottasch et al. ( 1984a ) and Lenzuni et al. (1989) were careful 
not to use fluxes from this band to derive the dust properties (or, 
rather, they used them onl y as up per limits). Some l ater st udies, 
however (Zhang & Kwok 1993| , Tajitsu & Tamura 1998) seem 
to have ignored this problem. 



Volk (1992) and Vil'koviskii & Efimov (1992) included the 



line e missio n in their model computations of the IRAS fluxes. 
Volk (1992), for example, showed that the IRAS color-color 



diagram was strongly affected by atomic line emission, prin- 
cipally in the 12 /im band. However, he explored only a very 
restricted parameter space, and it is interesting to see how the 
lines affect the IRAS fluxes in general during the planetary neb- 
ula evolution. 

Figure |l^-e shows the ratio Fi2(d+1)/Fi2(d) of the total 
emission in the 12 /im band (d+1 stands for dust + lines) and 
the emission due to dust only (d), as a function of S(lif3) for 
our model planetary nebulae. Figures ^-e, |^-e and ^-e are 
equivalent to figure [l|a-e for F25, Fqq and Fioo, respectively. 

We see that, as the nebula ages, the line contribution to the 
F12 emission increases. This is because, as the nebula expands, 
the grains become cooler and gradually emit at longer wave- 
lengths. With the gas chemical composition chosen in our mod- 
els, the lines that contribute most to the 12 /im band are [Ne II] 
12.8 /im, then [SIV] 10.5 /im, [ArlV] 9.0/tm, [ArV] 7.9 /tm 
and [Ne V] 14.3 /tm. The latter two lines become strong only if 
the central star is still luminous when it reaches a temperature 
above 150 000 K, which does not occurs for central stars with 
Af* > 0.62 Mq. Note that [Nelll] 15.5 /im, mentioned by Pot- 
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Fig. 1. The ratio Fi2(d-|-1)/Fi2(d) of total emission and pure dust emission in the 12 band vs. S(HP) for sequences of 
models computed at intervals of 500 yr until lOOOOyr. Each panel correspond to the variation of one parameter with respect to 
the reference model ( =0.60Mq, Mnob = O.2M0, Woxp = 20kms~^ and graphite grains with md/m^ = l .5 10^^). a: stellar 
mass M*: O.58M0 (circles), 0.60 (squares), 0.62Mq (triangles), 0.64 (diamonds), b: nebulai" mass M^oh'- 0.1 Mq 
(circles), 0.2 Mq (squares), 0.4Mq (triangles), c: expansion velocity Wcxp^ lOkms^^ (circles), 20kms^^ (squares), 40kms~^ 
(triangles), d: dust content md/m^: 7.5 10~^ (circles), 7.5 10~^ (squares), 7.5 10"'* (triangles), 7.5 10~^ (diamonds), e: dust 
material: graphite (circles), amorphous carbon (squares), circumstellar siUcates (triangles). 
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Fig. 2. Same as Fig. |l] but for the ratio F25(d+l)/i^25(d) 
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Fig. 3. Same as Fig. |l] but for the ratio F6o(d+l)/F6o(d) 
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Fig. 4. Same as Fig. |I] but for the ratio Fioo(d-|-l)/Fioo(d) 
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Fig. 5. The color temperature derived from the F25 / Fgo ratio 
as a function of the temperature derived by a mean square fit- 
ting of B\ (Td) to the observed emission in the three bands 
25, 60 and 100 ^m, giving equal weight to all three fluxes, for 
our sample of planetary nebulae. 



tasch et al. (1984a) does not contribute to the 12 /^m band. In 



the most extreme cases, lines may completely dominate the flux 
in the 12 /im band. Even for relatively young objects, the line 
contribution may still be important, depending on the charac- 
teristics of the object and its chemical composition. Therefore, 
unless a direct estimate of the infrared line strengths is avail- 
able, one should not use the observed flux in the 12 /im band 
(except as an upper limit) to infer average dust properties. A 
further problem with the 12 \iva band is that it is contaminated 
by the emission from the very small grains that experience tran- 
sient heating which could be especially important in the case of 
emission by policyclic aromatic hydrocarbons. 

It is interesting to note, in Fig.[l|d, that Fi2(d+1)/Fi2(d) is 
the smallest in the sequence of models with the lowest m^jm^. 
Intuitively, one might expect this ratio to be larger when the 
amount of dust in a nebula is smaller. What happens, in such 
a case, is that the total available energy to heat the grains is 
shared among a smaller number of particles, making each one 
hotter and the result is a larger dust emission at 12ptm. 

In comparison with the F12 emission, the line contribution 
to the ^25 emission is much smaller, but it can still dominate 
the dust emission by up to a factor about 3 in the later stages 
of evolution, as the lines contributing most are [OIV] 25.9 /im 
and [NeV] 24.3 ^m. 

The band with the smallest contribution of atomic lines in 
our models is the 60 ^m one (< 50 % in all cases), with the lines 
of [OIII] 51.8 ^m and [NIII] 57.3 /zm contributing mostly. 
This band is much wider than the 12 and 25 /im bands, and 
much of the dust emission in planetary nebulae goes through 
this band. 

The Fioo emission is not too much affected by line emis- 
sion ([OIII] 88.4 /im and [Nil] 121.9 /im being the contribu- 
tors) except at a very low dust-to-gas mass ratio (of about 10~^, 
see Fig.Hd). However, it has been reported (IRAS Explanatory 



Supplement 1988) that the 100 /zm band may be contaminated 



2 
2 .,-1 .,.-11 



by cirrus emission so that its use for planetary nebulae may be 
misleading unless a careful analysis of each case is made indi- 
vidually. 



log S(H/3) [erg cm s sr ] 

Fig. 6. The color temperature derived from the _F25 /^60 ratio 
as a function of 5(11/3) for our sample of planetary nebulae. 



5. The dust temperature 

Characteristic dust temperatures are generally obtained by fit- 
ting the observed spectrum with either B\{Td), where is 
the dust temperature, or a function which approximates the 
du st emis sivity: Bx {Td) or Bxj Td) ( e.g. Pottasch et 
al. [19844 Iyengar |1 984 Zhang & Kwok |l993| ). In the present 
study, we use A^^ B\{Td) because such a function better rep- 
resents the wavelength dependence of grain emissivity for the 
grains we have chosen. 

That such a quantity as a characteristic dust temperature 
can be defined in a nebula may be surprising. Indeed, there 
is a distribution of the temperatures of the dust grain particles 
inside the nebulae. Larger grains have smaller temperatures, 
since they are heated proportionally to their surface and radi- 
ate proportionally to their volume. Also, regions closer to the 
central star have higher dust temperatures because the stellar 
radiation field is less diluted. Finally, grains located in the ion- 
ized zone of the nebula are heated by a harder radiation field 
than those located in the neutral zone where only the stellar 
photons at energies below 13.6 eV and the Lya photons oper- 
ate. Nevertheless, our models show that, after integrating over 
all the grain particles in a given nebula, the resulting infrared 
spectrum is reasonably indicative of a characteristic dust tem- 
perature. The main parameter that affects the dust temperature 
during the evolution of a model nebula is the average distance 
of the dust particles from the central star, which varies by a 
factor 20 between 500 and 10 000 yr. 

To determine the characteristic dust temperature in the sim- 
plest way one can use the ratio of the infrared fluxes in two 
bands only, assuming that dust radiates like a single black- 
body with an emissivity proportional to A^^. This is a so- 
called color temperature. We have tested on our models that 
the color temperature derived from the F25(d)/^6o(d) ratio is 
very close to the temperature derived by a mean square fit- 
ting of A^^ B\{Td) to the pure dust emission in the 25, 60 and 
100/im bands, giving equal weight to all three fluxes (which 
confirms that the concept of a characteristic dust temperature 
is indeed a valid one). Similarly, these temperatures are also 
in agreement in the observed planetary nebulae, as shown in 
Fig.|^, where the contribution of lines (and ciiTus emission in 
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the band) adds some slight dispersion. Using a color 

dust temperature based on 25 and 60 fim fluxes allows us to de- 
rive characteristic temperatures in a consistent way for a large 
sample of planetary nebulae. From now on, Td will refer to the 
dust temperature determined in such a way. In principle, one 
should integrate A^^ B\{Td) over the IRAS bands but we have 
checked that doing this results in increasing Td by less than 
10%. Column 6 of Table 1 lists the values of Td for all the 
planetary nebulae of our sample. 

Figure ^ shows the values of Td as a function of the neb- 
ular average surface brightness, which has been computed us- 
ing the reddening corrected H/3 fluxes. There is a clear trend 
of decreasing Td with decreasing 5(H/3). This is reminiscent 
of the trend found by Pottasch et al. (1984ai and Lenzuni et 
al. (1989) of the dust temperature decreasing with increasing 
nebular radius, but the diagram has the advantage of being dis- 
tance independent. The dust temperatures in Fig. |] range be- 
tween 180K and 50 K. A few objects, namely M 2-23, M 1-48, 
He 2-34, PC 11, He 2-1 71, M 2-6 and H 1-36, have Td higher 
than 200 K and do not appear in the plot. 

Figure 0a-e shows the same as Fig. || for our sequences of 
model planetary nebulae, the values of Td being derived from 
the emission in the 25 and 60 /im bands including the lines. 
Roughly, the models show a trend of decreasing with de- 
creasing surface brightness. This is explained by a decreasing 
efficiency of the heating, mainly due to a dilution of the stel- 
lar radiation field inside the nebula as it expands. The models 
reproduce the observed diagram remarkably well, given their 
simplicity. Taking all the models together, one reproduces the 
observed range in Td, and the observed variation with S'(H/3). 
The increase in Td at the lowest S'(H/3) is due to an increasing 
contribution of the lines to F25, as seen in Fig. |[ 



6. The dust-to-gas mass ratio 



The formula used by Pottasch et al. ( |1984a| ) to derive dust-to- 
gas mass ratios in planetary nebulae depends on the assumed 
distance. We prefer using a formula that is distance indepen- 
dent. 

The infrared luminosity at the wavelength A is: 



AttcPFx 



ifr(a)pd(r,a)BA(Td(r,a))(4 7rr)^drda (1) 



-Ri, 



where: Fx is the measured flux from the source at the distance 
d; flmin and Omax are minimum and maximum dust grain sizes; 
K^^^{a) is the mass absorption coefficient (which, is almost 
independent of a at infrared wavelengths); pd{r, a) is the dust 
density. So, if Bx{Td) can be considered uniform in the nebula 
and independent of a, which, as we checked in our models, is 
not a bad approximation, then the dust mass is simply given by 



Fxd^ 



(2) 



The total luminosity in H/3 is approximately given by: 



47rd^FH/3 = / 4.110"22r<r"-^^nenH+47rr2dr (3) 

where the emissivity in H/3 is given by Eq. IV-25 of Pot- 
tasch ( 1984 ); Tc and ric are the electron temperature and den- 
sity, respectively; nH+ is the density of ionized hydrogen. If 
the gas is fully ionized and the density is constant the total gas 
mass is given by 



1.4mH47rd^FH/3 



4.1 10-22T 



-0.88„ 



(4) 



where mn is the mass of the hydrogen atom and the factor 1 .4 
comes from the classical assumption of a He/H ratio of 0. 1 . For 
simplicity, we will assume for the remaining of the paper that 

Te = 10*K. 

Therefore, using Eqs. and (|) the dust-to-gas mass ratio 
is given by 



Fx 



4.1 10-22T 



Fn0 Kf'BxiTd) 1.4mH47r 



(5) 



In principle, if the colour temperature is used, one can use 
equally well the fluxes in the 25 /im or the 60 fim band to es- 
timate m^/mg. However, since atomic lines may contribute to 
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the observed IRAS fluxes, we prefer to show results based on 
the 60 fim fluxes, since this band is less contaminated, as seen 
above (but the overall appearance of the diagrams based on the 
25 /im fluxes is very similar to the ones presented in this paper). 

To check the method before applying it to observed neb- 
ulae, we use it on our models, where we know the input dust 
properties and dust-to-gas mass ratio. Figure |^-e shows the 
values of m^/mg as derived from Eq. as a function of 
S{Y{j3) for our different series of models. In each case, we take 
the value of that corresponds exactly to the one used in 
the model: 145.32 cm^ g^^ for graphite, 74.38 for amorphous 
carbon and 53.45 for circumstellar silicates. 

We find that, for many of the models, the formula returns 
an estimated dust-to-gas mass ratio within ± 0.3 dex of the in- 
put one. Of course, the models that are optically thick return 
md/mg values larger than the input ones, since the H/3 flux 
corresponds to the ionized region only. Most of those optically 
thick models correspond to the vertical tracks seen in Fig. 8 
a-e, at the beginning of the planetary nebulae evolution. As 
seen in this figure, the total number of model points on these 
vertical tracks is small compared to the number of points in 
the later stages and they are located on the high surface bright- 
ness side of the diagram. Indeed, it can be shown that, whithin 
the theoretical framework adopted here, planetary nebulae with 
M* <0.62Mo are in an ionization bounded stage only dur- 
ing a short fraction of their observable lifetime, irrespective of 
the nebular mass and expansion velocity. For planetary nebule 
whose central stars have higher masses, however, the propor- 
tion of time spent in the ionization bounded regime is larger 
(see the figures in Stasinska et al. 1998), due to the fact that the 
stellar luminosity drop occurs early in the lifetime of the plan- 
etary nebula; the model with =O.64M0 in Fig. 8a, corre- 
sponds to such a case. There are however evidences that the to- 
tal number of planetary nebulae with central stars more massive 
than 0.62 M0 is not large, 20% at the very most (see Stasinska 
et al. 1997). We note also, in Fig.|^, and especially in panels a- 
c, that, for some of the very low surface-brightness models, the 
use of Eq. (^) underestimates the dust-to-gas mass ratio. These 
cases correspond to phases where the atomic line contribution 
to the 25 /im flux is important, as seen in Fig.^-e, so that the 
colour temperature of the dust is overestimated which induces 
an underestimation of the dust-to-gas mass ratio. In our series 
of models, this effect is really strong only in a few cases, mostly 
at the low surface brightness end. However, the effect will be 
larger at higher metallicities than considered in the models pre- 
sented here. 

It is interesting to compare the range in the derived md/mg 
values with the range of the input values of md/wg in our se- 
ries of models. Fig.^ shows that, if one excepts the optically 
thick cases, the method tends to narrow the md/mg range with 
respect to the input one at a given 5(H/3). This results from 
contamination of the IRAS bands by atomic lines, as can be un- 
derstood by considering Figs.^ and^. With larger abundances 
of the metals than in our models, the effect would be even more 
pronounced. This means that the method itself, when applied to 



log S(H/3) [erg cm s sr ] 



2 

2 ..-1 .,.-11 



Fig. 9. The dust-to-gas ratio derived from Eq. (|]) as a function 
of S(li(3) for our sample of planetary nebulae. 



real nebulae, does not artificially widen the distribution of dust- 
to-gas mass ratios with respect to S'(H/3). 

The sequence of models computed with circumstellar sil- 
icates systematically returns a smaller dust-to-gas mass ratio 
than the series computed with graphite or amorphous carbon 
grains (Fig.j^e). This is because in our determination of Td, we 
assumed that the dust emissivity varies like A^^, which is in- 
deed the case when comparing the 25 and 60 /im wavelenghts 
for graphite and amorphous carbon, while for circumstellar sil- 
icates the true dependence is rather like X~^. Therefore, we 
overestimate the dust temperature, and, consequently underes- 
timate the dust mass. This effect is not very strong, about a fac- 
tor 2 in m^/mg. However if we determine the dust-to-gas mass 
ratio of the model containing silicates using the value of Kqq^ 
for graphite, simulating the analysis of an object for which we 
do not know the true value of K^q^, the estimated dust-to-gas 
mass ratio goes down by an additional factor of about 2. 

In summary, combining all the uncertainties, the method 
to determine the dust-to-gas mass ratio should give estimates 
within, let us say, ±0.3 dex in the majority of planetary neb- 
ulae. As is known, it overestimates md/mg in the case of 
planetary nebulae that are incompletely ionized. Our models 
show that such cases, however, should represent a small frac- 
tion of observed nebulae, unless one is dealing with samples at 
high surface brightness. On the other hand, there is probably 
a tendency to underestimate the dust-to-gas mass ratio for low 
surface brightness objects, due to contamination of the IRAS 
fluxes by atomic lines, even when selecting the best IRAS 
bands, as we have done here. 

We now turn to the observational sample. Using Eq. (|]), 
where Fnp now represents the reddening corrected nebular flux 
in H/3 and taking for Kqq^ the value of 145.3 cm^g^^ cor- 
responding to graphite grains, we have computed the values 
m^/mg for our sample of planetary nebulae from the observed 
fluxes in the 60 /im band. These m^/mg are reported in col- 
umn 7 of Table 1. Figure ^ shows md/mg as a function of 
5(H/3). We find that in 70% of the cases the estimated values 
of md/mg lie between 10"^ and 10^'^, in 15% of them they 
are larger than 10~^ and in 15% of them they are smaller than 
lO^'^, the extremes being close to 10^^ and 10"^. Accounting 
for the fact that some planetary nebulae in the sample are ion- 
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Fig. 8. The values of ttla / rrig derived from Eq. (|]) for the models as a function of 5'(H/3). Same conventions as Fig. [l|. 



ization bounded, and that atomic line contribution may result in 
an underestimation of m^/nig in some low surface brightness 
objects, the true range in rrid/mg is probably of about a factor 
10. As a consequence, the dust-to-gas mass ratio is a primordial 
parameter for the interpretation of far infrared data of plane- 
tary nebulae. This was not always taken into account in former 
studies of statisti cal nat ure (Volk 1992 Zhang & Kwok 1993 , 
Tajitsu & Tamura[l998j). 

Clearly, in Fig. 9 there is no correlation between md/mg 
and S'(H/3) (the Pearson correlation coefficient is 0. 132 for 230 
points, with an associated probability of 0.05). Objects with 
very low or very high m^/m^ are found at any S'(H/3). The 
systematic errors in the determination of m^/mg that were 
discussed above would tend to overestimate m^/rrig at high 
5(11/3) and underestimate it at low S{i{(3). But, as we have 
argued above, important errors are not likely to affect a large 
fraction of the sample. We conclude that there is no evidence 
for a decrease in the dust-to-gas mass ratio as the planetary 
nebulae evolve. 

7. The total infrared flux and the infrared excess 

Several papers have studied the total infrare d flux a nd infrared 
excess in planetary nebulae. Pottasch et al. ( 1984a i not ed tha t 
younger planetary nebulae have higher IRE. Ratag et al. ( 199C ) 
compared a sample of about 100 planetary nebulae in the galac- 
tic bulge and a sample of about 100 bright nearby nebulae, and 
found that the latter showed, on average, lower IRE. They ad- 
vocated that one possibility could be for example that, because 
of selection effects, planetary nebulae from the bulge sample 
would be, on average, of higher surface brightness and thus, 
younger Another possibility suggested by Zijlstra et al. (1989) 
is that the higher IRE in the bulge planetary nebulae could be 
due to a higher dust content, related to a higher metallicity in 



the bulge. However, Ratag et al. (1990) rejected this option 
since the average metallicity of the bulge planetary nebulae is 
not that much different from the average metallicity of nearby 
disk planetary nebulae. They then argued that the bulge plan- 
etary nebulae, being descendent of a lower mass stellar pop- 
ulation, would have lower mass nuclei. These would spend a 
longer time at lower effective temperature, inducing higher IRE 
due to the larger contribution of stellar non ionizing photons to 
the dust heating. Finally, they explored the possibility that the 



-9 

log F|„ 



(25,60) [Jy] 



Fig. 10. The total infrared flux estimated from the IRAS fluxes 
at 12, 25, 60 and 100 /im as a function of the total infrared flux 
estimated from the IRAS fluxes at 25 and 60 /im only 



high IRE in the galactic bulge planetary nebulae could be due to 
additional dust heating by the interstellar radiation field. They 
estimated that, for a planetary nebula located in the galactic 
bulge, an energy of 100 to 400 would be provided by the 
interstellar radiation field. 

Our models are able to provide some insight into these 
questions, but let us first discuss the empirical determination 
of the total infrared flux. 

The total radiation that is emitted by dust in a nebula can 
be estimated from the IRAS fluxes basically by two methods. 
One is simply to compute the flux under the curve defined by 
the observed infrared fluxes, as done by Pottasch et al. (1984a). 
The other is to fit the observed energy distribution by a func- 
tion proportional to A^^ B\{Td)- The simplest version of this 
method is to compute the blackbody flux that corresponds to 
the ratios of the observed fluxes at 25 and 60 /im. Beside the 
advantage of being computationally simple, this method allows 
one to consider a large number of planetary nebulae, since it 
requires fluxes in only two bands, and it is the least affected 
by atomic line and cirrus emission. It is interesting to compare 
the results of these two methods. This is done in Fig. where 
we plot the total IRAS flux determined by the first method us- 
ing the fluxes at 12, 25, 60 and 100 /im as a function of the 
total IRAS flux derived by the simplest version of the second 
method. In building this diagram, we considered only those ob- 
jects with Q = 3 in the four IRAS bands, so that the total num- 
ber of objects represented is only 107. The observed correlation 
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Fig. 11. The values of IRE derived from Eq. (|6|) as a function 
of the angular distance to the Galactic center for our sample of 
planetary nebulae. 
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Fig. 12. The values of IRE derived from Eq. (|6|) as a function 
of SiHP) for our sample of planetary nebulae. 



is very good over the whole range of fluxes. From this, we feel 
confident that we can indeed use only the 25 and 60 /im fluxes 
to estimate the total infrared flux. The infrared excess is then 
determined from the total infrared flux and the total reddening 
corrected H/3 flux using 



IRE = 



Ft 



IR 



23.3 Fh/3 



(6) 



where 23.3 is the ratio of Lya to H/3 emission (Table VIII-2 of 



Pottasch |1984 i. 

Figure |l ifshows the IRE, determined in such a way, as a 
function of the angular distance to the Galactic center for our 
sample of planetary nebulae. The sample is much larger than 
considered in previous studies since it contains more than 500 
objects. There is indeed a tendency for the average infrared ex- 
cess to decrease with increasing angular distance to the Galac- 
tic center. However, the variation is not exactly the one that 
would be expected if the interstellar radiation field alone were 



to produce this effect, as suggested by Ratag et al. (1990) and 



Amnuel (1994). Indeed, the variation is very smooth and shows 



no marked peak in the direction of the bulge. We also note the 
existence of objects (about 10% of the whole sample) with IRE 
> 10, and those are distributed uniformly in Galactic longitude. 
This indicates that these high IRE are related to an intrinsic 
property of the nebulae and not to their location in the Galaxy. 



Figure [T^ shows the IRE as a function of Silip). In this 
figure, which contains over 470 points (the angular diameter is 
required in addition to the H/3 flux), one can see that planetary 
nebulae of high or low surface brightness can have equally high 
IRE. There is no evidence that the IRE is higher for the younger 
nebulae. 

Figure 13 i-e represents the same diagram as Figure u2 for 



our sequences of models. It shows that, for the earliest evo- 
lutionary stages, when the central star does not emit a large 
proportion of ionizing photons, the infrared excess can indeed 
reach very high values, as was argued in several previous stud- 
ies. But, according to our models, such optically thick objects 
should be rare. On the other hand, large values of the IRE are 
easily produced when m^/mg is high (Fig. pj|d). Since we have 
shown that large dust-to-gas ratios occur for planetary nebulae 
of all ages, the same is true for large values of the IRE. 

It has been suggested by Zijlstra et al. ( 1989| ) that the in- 
frared excess may be used as a measure of the temperature 
of the central stars, Tcff, for stars cooler than 40000K (for 
larger values of Toff, the proposed calibration becomes mul- 
tivalued). Figure [ija-e shows the IRE as a function of Toff 
in our sequences of models. From this figure, we see that the 
method is less reliable than it appears from Fig. Bl of Zijlstra 
et al. (1989). Even for the coolest stars, for a given Tcs, the 



IRE depends somewhat on the mass of the central star, on the 
nebular mass and expansion velocity and, most importantly, on 
the dust-to-gas mass ratio. As a result, the derivation of Teff 
using the infrared excess is not accurate. It will be especially 
misleading if the dust-to-gas mass ratio in a given nebula does 
not al low the simplifying assumptions made by Ziljstra et al. 
(|1989|) to be valid. 



8. The grain size 



Lenzuni et al. (1989) estimated the typical size of the dust 
grains in their sample from the thermal balance equation, un- 
der the assumption that the absorption efficiency Quv was 1 
in the UV and that they could estimate the total luminosity of 
the central star. They found that the average size of the grains, 
a, was decreasing from about 3. 10^^ cm to 3. 10^^ cm as the 
radius of the nebulae increased from 3. 10^^ to 1. lO^'* cm. This 
was an appealing conclusion, which, combined with their find- 
ing that the total dust-to-gas mass ratio was decreasing with 
time, led to some interesting speculations about the destruction 
process of dust grains in planetary nebulae. 

Actually, it is not quite true that Quv is constant, it varies 
strongly with wavelength for grains smaller than about 10^^ - 
5. 10~^ cm depending on the dust material. 



The method of Lenzuni et al. (1989) can easily be tested on 



our models, which all have the same distribution of grain sizes 
during the whole nebular evolution. The formula they used can 
be rewritten as: 



3 A ^ Q\]w 



bulk 



(7) 
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Fig. 13. The values of IRE derived from Eq. (||) for the models as a function of 5(H/3). Same conventions as Fig. |l|. 
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Fig. 14. The values of IRE derived from Eq. for the models as a function of Toff. Same conventions as Fig. |l| 




where L» is the stellar luminosity and Pj""^ the density of 
the grain material. In the case of our models, we know the 
total stellar luminosity, so the formula can be applied easily. 
Fig. |l5^-e shows the resulting value of a obtained using Eq. (0) 
at A = 60 fim for our sequences of models. We see that the 
method to derive a is strongly biased, and gives decreasing 
values of a as 5(H/3) decreases. The effect is quite large: the 
estimated a decrease by about 2 dex as S{ll(3) decreases by 
about 4 dex which, in our models, correspond roughly to an in- 
crease of the radius by about 1 dex. Therefore, the bias in the 
method is of the same order of magnitude as the effect found 
by Lenzuni et al. ( |1989| ). 



We h ave not questioned here the way in which Lenzuni et 
al. ( 1989 ) determ ined L ^,. They used an analytical model of HII 
regions (Panagia 1974 ) to do this. This model involves many 



assumptions that are not necessarily fulfilled and might pro- 
duce an additional bias in one way or another. 

Since, in our models, the value of a is 3.53 10^^ cm, we 
see from Fig. |l5| that the method of Lenzuni et al. rather tends 
to underestimate a. If we consider that the highest values of 
a found by Lenzuni et al. are ~ 3. 10~^ cm, it might be that 
indeed such large grains are found in planetary nebulae. Such 
a conclusion, however, might not be robust, since we have not 
tested how it would vary if we assumed different size distri- 
butions or forms of the grains. Interestingly, we note that, in 
the Red Rectangle nebula. Jura ( 1997 ) has argued, on other 
grounds, that grains have a large size. 

The main conclusion from testing Eq. (^ on our models is 
that, clearly, one cannot derive the value of a in such a way. One 
could, perhaps estimate this parameter from a detailed model 
fitting of individual, well observed, objects, but this is not cer- 



G. Stasinska & R. Szczerba: The dust content of planetary nebulae 



11 



tain and has never been attempted so far to our knowledge. For 
the time being, we consider that there is no evidence that the 
typical grain size of dust varies during the evolution of plane- 
tary nebulae. 

9. Discussion 

We have reconsidered the dust content of planetary nebulae by 
performing a new statistical study of broad band IRAS data. 
The main difference with previous studies is that our method of 
analysis was tested on a wide grid of photoionization models 
of evolving dusty planetary nebulae. Furthermore, our results 
are based on distance independent diagrams. We are aware that 
our analysis, although more elaborated than previous studies of 
statistical nature, still suffers from uncertainties. Detailed mod- 
elling of individual objects should, in principle, provide better 
estimates of the dust content of planetary nebulae. Such an ap- 
proach is however not feasible on a large sample of objects. 

Our main result is that the range in dust-to-gas mass ratios 
among planetary nebulae is large, about a factor 10, but, in con- 
trast with previous authors, we find no evidence for a change 
in the dust-to-gas mass ratio in the course of planetary nebu- 
lae evolution. Neither do we find evidence for a change in the 
average grain size. 

To our knowledge, the problem of dust survival in plane- 
tary nebulae is still poorly understood from a theoretical point 
of view. The main processes of dust destruction have been de- 



scribed by Draine & Salpeter (1979 and references therein) or 



Draine (1990) but the estimated efficiencies are uncertain and 



need confirmation by observation. After the works of Pottasch 



et al. (1984a) and Lenzuni et al. (1989) the idea that dust is 



gradually destroyed - by the energetic photons emitted by the 
central star or by other mechanisms - made its way among 
astronomers (e.g. F urton & Witt 1992 , Dinerstein et al. 1995, 
Corradi & Schwarz 1995 ). Our analysis does not support such 
a view. In addition, there is now observational evidence that 
dust destruction in planetary nebulae i s prob ably not as strong 
as thought before. Namely, Volk et al. ( 1997 ) have reported that 
their attempt to search for lines of Ca in planetary nebulae has 
failed, indicating that Ca is probably locked in grains (see the 
classical work of Field 1974 which discusses the problem of the 
gas-phase element depletions) and the UV radiation does not 
seem to be so efficient for dust destruction. It has been shown 



by the theoretical calculations of Okorokov et al. (1985) and 
Marten et al. (1993) that radiation pressure on dust particles 
accelerates particles to higher velocities and gradually expells 
them from the nebula. It seems however, that the conclusion 



of Okorokov et al. (1985) that dust is removed from planetary 



nebulae almost entirely is too far going. They did not consider 
the process of charging of the dust particles which probably 
freezes them into the ionized gas. It is worth of noting, how- 
ever, that the spread of the points representing the dust-to-gas 
mass ratio is quite wide and shows that planetary with differ- 
ent dust contents exist. The reason could be at least twofold: 
planetary nebulae could be created with different dust-to-gas 
mass ratios on the AGB; or the time spent by the star during 



the pre-planetary phase of evolution, when the process of dust 
acceleration and removing is efficient, was different (shorter 
for planetary nebulae with presently higher dust-to-gas mass 
ratio). If the second scenario applies, one can speculate that, 
at least from a statistical point of view, planetary nebulae with 
smaller dust-to-gas ratio might have less massive central stars. 
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